Monte Carlo EM for Generalized Linear Mixed Models using Randomized Spherical Radial Integration

The expectation-maximization algorithm has been advocated recently by a number of authors for fitting generalized linear mixed models. Since the E-step typically involves analytically intractable integrals, one approach is to approximate them by Monte Carlo methods. However, in practice, the Monte Carlo sample sizes required for convergence are often prohibitive. In this paper we show how randomized spherical-radial integration (Genz and Monahan, 1997) can be implemented in such cases, and can dramatically reduce the computational burden of implementing EM. After a standardizing transformation, a change to polar coordinates results in a double integral consisting of a one dimensional integral on the real line and a multivariate integral on the surface of a unit sphere. Randomized quadratures are used to approximate both of them. An attractive feature of the randomized spherical-radial rule is that its implementation only involves generating from standard probability distributions. The resulting approximation at the E-step has the form of a fixed effects generalized linear model likelihood and so a standard iteratively reweighted least squares procedure may be utilized for the M-step. We illustrate the method by fitting models to two well-known data sets, and compare our results with those of other authors.